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study the dynamics of thermonuclear flames propagating in fuel stirred by stochastic forcing. 
ii^The fuel consists of carbon and oxygen in a state which is encountered in white dwarfs close to 
^he Chandrasekhar limit. The level set method is applied to represent the flame fronts numerically. 
^VThe computational domain for the numerical simulations is cubic, and periodic boundary condi- 
C _ } tions are imposed. The goal is the development of a suitable flame speed model for the small-scale 
QQiynamics of turbulent deflagration in thermonuclear supernovae. Because the burning process in a 
f~*y upernova explosion is transient and spatially inhomogeneous, the localised determination of subgrid 
i^~£cale closure parameters is essential. We formulate a semi-localised model based on the dynamical 
.^ equation for the subgrid scale turbulence energy fcsgs- The turbulent flame speed st is of the order 
■^^y^2A;sgs. In particular, the subgrid scale model features a dynamic procedure for the calculation of 
r^ he turbulent energy transfer from resolved toward subgrid scales, which has been successfully ap- 
Qp ^ied to combustion problems in engineering. The options of either including or suppressing inverse 
I energy transfer in the turbulence production term are compared. In combination with the piece-wise 
Qparabolic method for the hydrodynamics, our results favour the latter option. Moreover, different 
^Hchoices for the constant of proportionality in the asymptotic flame speed relation, st tx ■^2fcags, are 
'^investigated. 

^Keywords: Combustion, thermonuclear, turbulence, large eddy simulation, level set method 
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Introduction 



A certain kind of stellar explosion, known as type la supernovae among 
astronomers, is currently explained by the thermonuclear explosion of an 
electron-degenerate stellar remnant Such an object, which is called a white 
dwarf, emanates from the burn-out of stars comparable in mass to our Sun and 
is mainly composed of carbon and oxygen. If the white dwarf has a compan- 
ion star in close orbit, it can grow by accreting material from the companion. 
Under certain conditions, the white dwarf's mass will steadily increase and 
finally approach the Chandrasekhar limit, which is the maximal mass that can 
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be supported by the degeneracy pressure of electrons As the temperature 
and density are increasing, thermonuclear burning of carbon and oxygen grad- 
ually sets in. Close to the Chandrasekhar mass, the conditions in the core of 
the white dwarf eventually pass a critical threshold [HJ. At this point, the rate 
of thermonuclear reactions rises rapidly, and a runaway is initiated, which in- 
cinerates and disrupts the whole star within a few seconds. The total energy 
release is of the order 10^^ erg j^. 

The thermonuclear combustion of degenerate carbon and oxygen of density 
in the range ~ 10^.. 10^ gem proceeds in the form of a deflagration [21 
10]. Since the nuclear ash produced by the burning process has less specific 
weight than the surrounding unprocessed material, it becomes Rayleigh- Taylor 
unstable. Since turbulence is subsequently produced, the flames get corrugated 
and folded In consequence, there is a positive feedback mechanism of 

turbulence enhancing the burning and, according to state-of-the-art numerical 
simulations, eventually results in an explosion [9|lin|[TT]. Although it cannot 
be ruled out that a transition from the deflagration to a detonation might 
set in at some stage fT2|[T3]. turbulent deflagration plays a crucial role in the 
theoretical modelling of thermonuclear supernovae in any case. 

The subject of this article is the dynamics of flame fronts on length scales 
much smaller than the size of a Chandrasekhar-mass white dwarf. To that 
end, an artificial scenario was set up. A turbulent flow is produced by means 
of stochastic stirring in a cubic domain subject to periodic boundary condi- 
tions |141ll5j . Thermonuclear burning is ignited in small spherical regions and 
subsequently evolved by means of the level set method 1 6 17^ . The complicated 
network of thermonuclear reactions encountered in a type la supernova is sub- 
stituted by the effective fusion of equal mass fractions of ^^C and ^^O to ^^Ni 
and ^He as representative reaction ^H] . The equation of state is dominated by 
the degenerate gas of relativistic free electrons. Thus, the approximate relation 
P oc yO^/^ applies, while the temperature has virtually no influence on the pres- 
sure. This is actually the reason for the runaway, because the negative feedback 
between heating and expansion in non-degenerate matter is absent. The exact 
equation of state has no analytic solution and must be integrated numerically. 
Moreover, contributions from nuclei, photons and pair electron-positron pair 
creation at temperatures of the order 10^'' K are taken into account (section 
3.2 of 19 ). The fluid dynamics is treated by means of the piece-wise parabolic 
method (PPM) within the framework of the Euler equations j2U| . 

In the corrugated flamelet regime of combustion, the flame propagation is 
affected by turbulence on length scales ranging from the Gibson length up to 
the integral length scale [25 . In general, only the the largest length scales can 
be resolved in a numerical simulation. In order to account for the wrinkling of 
the flame surface on length scales smaller than the cutoff scale of a simulation, 
an effective propagation speed, the so-called turbulent flame speed, must be 
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calculated. This involves a subgrid scale (SGS) model for the local budget of 
turbulence energy contained in numerically unresolved modes. In this article, 
we present a numerical study in which different variants of the SGS turbulence 
energy model are compared (section 4.3 in |12])- In particular, we adopted a 
dynamical procedure for the computation of SGS closure parameters. This 
procedure was proposed by Kim and Menon for the application in LES of gas 
turbine combustor flows |23j . 

2 The physics of turbulent thermonuclear deflagration 

The mechanism of deflagration is based on thermal conduction, as opposed 
to a detonation, which proceeds via shock compression. Unburned material 
(fuel) is heated in the vicinity of the reaction zone and thereby gets ignited. 
Once heat generation is balanced by diffusion, the burning zone is propagat- 
ing at a steady subsonic speed, and pressure equilibrium is maintained across 
the reaction zone. Basically, this characterises what is commonly known as a 
flame. For chemical combustion, a distinction is made between premixed and 
diffusive flames. Thermonuclear flames are trivially premixed, because no ad- 
ditional agent, like oxygen in most chemical burning processes, is required. 
The local propagation speed of the flame, which is solely determined by mi- 
croscopic properties, is called the laminar burning speed. The notion of a flame 
applies, if fluid motions do not significantly disturb the burning process within 
the reaction zone, i. e., the characteristic time scale of burning is much smaller 
than the kinetic time scale of velocity fluctuations on length scales compara- 
ble to the flame thickness 5-p. Equivalently, 5-p *C /g) where the Gibson length 
scale Iq is the smallest length scale on which the burning process is affected 
by fluid motion |2J . The condition 8y <C thus specifies the flamelet regime 
of combustion, which is reviewed in this section. For the thermonuclear com- 
bustion in C+0 white dwarfs, it appears that the flamelet description is valid 
for /O > 3 • lO^gcm"^ ^Ij. In thermonuclear supernovae, most of the burning 
takes place at significantly higher densities. 

2.1 Laminar burning 

The width of the reaction zone, (5f, is determined by the equilibrium between 
energy generation due to nuclear reactions and the rate of diffusion caused 
by thermal conduction (§ 128 in |2S])- The balance between these processes 
can be expressed in terms of their characteristic time scales, Tburn and Tcond- 
The former is given by Tburn ~ P^nuc/B, where e^uc is the energy generated by 
the fusion of a unit mass of nuclear fuel, and B is the rate of energy release 
per unit volume. The time scale of conduction, on the other hand, can be 
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expressed as Tcond ~ S^/lcC, where le is the mean free path of the electrons, 
which contribute the major part of the thermal conductivity, and c is the 
speed of sound. Setting these two time scales equal, one finds that the flame 
thickness 5f is approximately given by 



t'F ~ Y ^ • (1) 

Defining the laminar flame speed by siam = ^p/T'burnj we have 



/ ^^^^ 

•Slam ~ -\/ • 

V /Offnuc 

The specific energy release for the fusion of "'^^C and "'^^O to ^^Ni is £"nuc ^ 
7 ■ lO^^ergg"^ ^H]. The flame speed siam for the thermonuclear combustion 
of degenerate carbon and oxygen was computed numerically for a wide range 
of mass densities and nuclear compositions by Timmes and Woosley [Oj. For 
example, siam ~ 3.6- 10® cms~^ and dp ~ 2.9- 10"'^ cm for equal mass fractions 
of carbon and oxygen at a density 10^ gcm~^. 



2.2 Turbulent burning 

So far, we have only been concerned with the microphysics of thermonuclear 
deflagration. Let us now consider the combustion of C+0 fuel in a state of 
turbulent motion. For brevity, we shall assume the case of steady isotropic 
turbulence, i.e. a statistically self-similar hierarchy of vortices or eddies. Each 
vortex of size I has a characteristic velocity v'(l) and an associated turn-over 
time Teddy (0 = l/v'{l). If v'{l) is small compared to the laminar flame speed 
■Slam) then the flame front will propagate through a region of size I in a time 
much faster than the turn-over time reddy(0- Hence, the turbulent flow appears 
to be more or less "frozen" with respect to the burning process on these scales. 
For v'{l) ^ Slam; on the other hand, the front is significantly distorted while 
it is crossing a vortex of diameter /. Hence, there is a threshold scale on which 
burning decouples from turbulence. This is the Gibson length Iq, which is 
defined by [21| 

v'{Ig) = Slam- (3) 

At length scales I ^ Iq, turbulence corrugates the flame and thereby increases 
its surface area. Consequently, turbulence enhances the burning process and 
the release of heat is growing. This can be accounted for by introducing a 
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turbulent propagation speed St{l). In other words, averaging over regions of 
size I, the flame front propagates with an effective speed st{l) greater than the 
laminar burning speed siam, which specifies the local speed of any portion of 
the flame smaller than ^g- 

For I ^ Iq, st{l) becomes asymptotically independent of siam- The funda- 
mental hypothesis applied in this article is that st{l) is then given by the 
magnitude of the turbulent velocity fluctuations v'{l). In the framework of the 
phenomenological Kolmogorov theory of isotropic turbulence this veloc- 
ity obeys the scaling law v'{l) oc l^^'^ in the inertial subrange, i.e. in the range 
of length scales which are neither affected by the viscosity of the fluid nor 
large-scale energy injection. Therefore, 

st{l) ~ v'{l) oc /V3. (4) 

The relation St{l) ~ was first proposed by Damkohler, who studied Bun- 
sen cones in the laboratory |^7j. Further validation of this conjecture came 
from numerical studies |2S|- A motivation based on a theoretical analysis in 
the framework of the level set prescription was given by Peters |29j . 

3 The numerical modelling of turbulent flame propagation 

In early studies of thermonuclear deflagration [3101, a reactive-diffusive flame 
model with artifical diffusion and reaction rates was applied. In this approach, 
the thickness of the flame is artifically increased over several grid cells and 
the propagation speed is adjusted to a prescribed value. On the other hand, 
the level set method proposed by Osher and Sethian [THIIHT] is a front tracking 
method which describes the interface separating ash from fuel as a genuine 
discontinuity. The interface is numerically represented by the set of all points 
for which a suitably chosen distance function vanishes, i.e. the zero level set. 
This is a sensible approximation if the physical flame thickness is very small 
compared to the Gibson scale. For the simulation of thermonuclear combustion 
in type la supernovae, the level set method was implemented by Reinecke 

[mill]. 

3.1 The level set method 

Let G{x,t) be a signed distance function with the property |VG| = 1. The 
absolute value |G(a;,i)| is equal to the minimal distance of the point x from 
the flame front at time t. The front itself is given by the constraint G{x, t) = 0, 
i.e. it is represented by the zero level set T{t) = {x\G{x, t) = 0}. With the sign 
convention G{x,t) > in regions containing burned material, the unit normal 
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vector pointing towards unburned material is given by n = — VG/|VG|. The 
time evolution of the front T{t) is implicitly determined by the total time 
derivative of G{xrit),t) = 0. For a certain point at the front, xr{t) G r(t), we 
have 



The speed function is given by the sum of two contributions. Firstly, the 
advection speed normal to the flame front, Vy^ ■ n, where i»u is the velocity of 
the fuel immediately ahead of the front in an Eulerian frame of reference. And, 
secondly, the intrinsic propagation speed s of the flame front relative the fuel. 

The local equation can be formulated globally as well, without con- 
straining the position x to the flame surface. Substituting the definition of the 
normal vector n and expressing the speed function in the form v^^ - n + s, the 
evolution equation for the level set function at any point in space becomes 



The advection part on the right-hand side can be treated with a finite-volume 
scheme, for instance, the PPM. The intrinsic front propagation is usually calcu- 
lated by means of an entropy-satisfying upwind scheme. In general, non-planar 
fronts will develop sharp corners and the corresponding level set must be a 
weak solution: Information about the initial conditions is lost, once a cusp 
has formed, and the subsequent evolution is irreversible. The corresponding 
entropy condition can be formulated in the following way: Once a certain fluid 
element is burned, it remains burned thereafter. In fact, this implies the equiv- 
alent Huyghen's principle in optics for the propagation of the front over an 
infinitesimal interval of time (section 5 in Finally, in order to preserve 

the property \'VG\ = 1, the updated distance function has to be corrected 
after each time step. In the implementation of Reinecke, this is achieved by 
means of re-initialisation |17j . 

For the complete implementation of the level set technique, both the burned 
and the unburned state in an intersected numerical cell must be reconstructed 
from the jump conditions across the front. Assuming that there is a volume 
fraction of unburned material a, conservation of momentum imposes the con- 
straint 





dG{x,t) 
Ft 



[v^{x,t) + s{x,t)n{x,t)]\VG{x,t)\. 



(6) 



pv = apnV^ + (1 - a)PhVh 



(7) 



given the finite- volume averages p and v. The volume fraction a can be cal- 
culated by linearly interpolating the discrete numerical values of the distance 
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function G. Supplementing the momentum equation with the Rayleigh cri- 
terion, the Hugoniot jump condition and the continuity constraint for the 
tangential velocity components, a non-linear system of equations is obtained. 
The solution yields v^, and the corresponding state variables (cf. [32|). 
This procedure of in- cell reconstruction was indeed successfully implemented 
for chemical combustion problems 33 . However, the deviation of the interpo- 
lated front element from the exact smooth solution can introduce significant 
errors in the reconstructed states. In particular, it is sometimes impossible 
to reconstruct physically sound states for degenerate matter, because of the 
stiffness of the equation of state. Moreover, one faces topological ambiguities 
for certain configurations. A pragmatic method is to average over all possible 
values, whenever one of these rare cases is encountered fI7|. Although Ropke et 
al. have recently succeeded with the implementation of in-cell reconstruction 
for the problem of thermonuclear flame propagation in two dimensions |32| . 
generalising the algorithm to three dimensions would be much more challeng- 
ing. 

The difficulties outlined above are avoided with the so-called passive imple- 
mentation, where the difference between burned and unburned states is ne- 
glected and the advection speed is set equal to v n. The discrete values of the 
velocity and state variables are then interpreted as cell-centred averages. This 
is a fair approximation in the limit of moderate density jumps between fuel 
and ash. A caveat of using the passive implementation for simulations of burn- 
ing at low density is the generation of numerical artifacts. Fortunately, these 
problems are mainly encountered for densities significantly less than about 
10*^ gcm~^. Apart from the systematic errors introduced by the averaged den- 
sity and advection velocity, the burning zone is not strictly represented by the 
zero level set. Actually, there is a mixed phase between the regions contain- 
ing pure ash and fuel, respectively. The width of the diffusive smearing of the 
flame is typically a few cells, which is still less than for the reaction-diffusion 
method. It was demonstrated by numerous applications in simple test prob- 
lems as well as large-scale simulations of thermonuclear supernovae that the 
passive implementation gives a satisfactory representation of the flame fronts 
at high density and is robust even in three dimensions |34|lin). For this rea- 
son, we used the passive implementation for the simulations presented in this 
article. 

3.2 The turbulent flame speed relation 

On length scales larger than the Gibson length Iq, flames are predominantly 
shaped by turbulence. If Iq is significantly smaller than the numerical res- 
olution A, the computed flame front appears inevitably smoother than its 
physical counterpart. Consequently, the predicted burning rate would be un- 
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derestimated, if just the laminar burning speed was substituted for the intrinsic 
propagation speed s in equation ©. This is where the notion of the turbulent 
flame speed comes in. We propose that St(A) is given by the local magnitude 
of unresolved turbulence velocity fluctuations v'{A) [2], albeit the turbulent 
flame speed, in a strict sense, is an ensemble average. Since the propagation 
speed cannot be less than the laminar flame speed, the simplest relation with 
correct asymptotic behaviour is 



st(A) = max(siam, \/2C^A~) = max(siam, yC^^sgs)- (8) 

The quantity /Cggs = l^sgs the subgrid scale turbulence energy and gsgs ~ 
v'{A) the corresponding speed. An exact definition will be given in next sec- 
tion. For brevity, it is understood that st denotes the turbulent flame speed 
at the numerical cutoff A in the following. 

For A ^ Iq, we have the asymptotic relation st ~ VC^Qsgs in the limit 
of fully developed turbulence. Consequently, the turbulent flame speed be- 
comes independent of the laminar flame speed, and the parameter y/Ct de- 
termines the asymptotic scaling of the turbulent flame speed. However, it is 
not quite clear whether the constant of proportionality in the relation between 
St (A) and gsgs is just unity or a different value. Empirically, it appears that 
St(A) = 2v'(A), where v'{A) = Qsgs/V^ EHI- Thus, Ct = 4/3 in agreement 
with a constant of proportionality close to unity. Another shortcoming of the 
maximum relation Q is that it gives a good approximation to the turbulent 
flame speed in the laminar and the fully turbulent regime, respectively, but 
not for the transition in between. If gsgs ~ siamj the relation between turbulent 
flame speed and turbulence velocity might very well be different. For exam- 
ple, Im et al. mention a quadratic dependence on the turbulence velocity in 
the case of weak turbulence I35j . On the other hand, Ropke et al. report a 
linear relation even for turbulence velocities which are only marginally larger 
than the laminar flame speed in two-dimensional numerical simulations j^fij . 
However, this result is possibly unsubstantial for the three-dimensional case. 
Apart from that, the transition from laminar to turbulent burning progresses 
rather quickly, and a correct description in the intermediate phase is therefore 
not overly important. 

A different turbulent flame speed model was motivated theoretically by 
Pocheau [37; : 

(\ n-| l/n 
^ . (9) 

■5lam / 

In the scale- invariant regime, with q^gs ^ siam, the asymptotic form st ~ 



■St 
■^lam 
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Ct (Zsgs is obtained. Kim et al., chose n = 2 and Ct ~ 20/3 for LES of gas 
turbine combustor flows '7T . This value was inferred from several laboratory 
experiments with hydrocarbon/air flames. However, the data points cover val- 
ues of the turbulent flame speed of the same order as magnitude as the laminar 
burning speed only. For this reason, one must be careful with any extrapolation 
to the fully turbulent regime, in which st ^ siam- If Qsgs ^ •sianu Taylor series 
expansion of the right-hand side of equation © yields the Calvin- Williams 
relation, 



which is consistent with the numerical results of Im et al. pS] . 

Apart from calculating the turbulent flame speed, secondary SGS effects can 
be included in the dynamical equation © for the level set function j23|l29j . 
This was indeed numerically investigated by Im et al. In particular, they 
suggested a procedure for the computation of Ct in the fashion of the localised 
closure for the production parameter C^, which will be discussed in section 
Whether this is advisable in combination with the passive implementation, 
where numerical artifacts in the shape of the resolved level set might produce 
significant spurious contributions, is questionable. For this reason, it has not 
been attempted. In addition, there is a SGS transport term for the level set, 

which is of the form dk{{vkG)cs — v^G) and effectively introduces diffusion 
of the level set due to SGS turbulence. However, Kim et al. argued that the 
contributions arising thereof are not particularly important and, in fact, cannot 
be determined within the available framework of SGS modelling |23j . 



4 The subgrid scale model 

For the determination of the turbulent flame speed according to equation (jSJ 
or (jUJ, the kinetic energy k^gs of unresolved vortices has to be computed. A 
dynamical equation for /cggs is obtained through decomposition of the conser- 
vation law for kinetic energy. The procedure of decomposing is conceptually 
based on the notion of filtered quantities. In general, a filter is a convolution 
operator, which smoothes out fluctuations on spatial scales smaller than the 
characteristic length of the fllter. If a certain numerical solution of the hydro- 
dynamical equations is computed by means of a finite- volume scheme, say, the 
PPM, then one can associate this solution with the smoothed velocity field 
v{x,t), which is obtained by mass- weighted or Favre filtering of the exact 




(10) 
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realisation of the flow: 



v{x, t) 



p{x,t) 



(11) 



where p{x,t) = {p{x,t))eff is the smoothed mass density. The underlying 
hypothesis is that, if the physical flow 'v{x, t) were known, there exists a filter 
( )eff with suitable properties such that the Favre-filtered velocity field v{x,t) 
would reproduce the numerically computed velocity field. Corresponding to 
the smoothed and the fluctuating components of v{x,t), respectively, one can 
distinguish the resolved part, fcres = ll"*^!^) and the subgrid scale part kggs 
of the specific kinetic energy. In the following, a formal decomposition of the 
kinetic energy is devised and the dynamical equation for kggs is formulated. 
The non-linearity of the conservation laws necessitates closure relations for 
several terms in the decomposed equations. SGS closures and the calculation 
of associated parameters are discussed in the remainder of this section. 

4.1 The subgrid scale turbulence energy model 

In Germano's consistent decomposition, the SGS turbulence energy is simply 
defined by the difference between smoothed and resolved kinetic energy |22| 
155] . This decomposition is equivalent to setting /cggs = —^tu, where the SGS 
turbulence stress tensor r^fc is defined by 



We prefer the Germano decomposition, because it yields the conceptually most 
transparent definition of the SGS turbulence energy and avoids formal difficul- 
ties associated with SGS closures. Deviations of the turbulent flame speed in 
alternative decompositions correspond to higher-order terms which result from 
secondary filtering of filtered quantities [22]. These contributions are likely to 
be insignificant within the intrinsic inaccuracy of the flame speed model 

The SGS turbulence stress tensor Tik enhances the viscous dissipation given 
by cTjfc in the equation of motion for the filtered velocity field: 




(12) 



Hence, the SGS turbulence energy is given by 




(13) 




(14) 
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Here it is assumed that the specific stirring force f- injects energy on length 
scales which are large compared to the numerical resolution A. The operator 
D/Dt is the Lagrangian derivative, 



The dynamical equations for kinetic energy in the Germano decomposition 
read 



P j-j^ "-res — Vi 



dP (s)^ d 



(16) 



Pj^^^sgs -^sgs — ^sgs p('^sgs ~l~ £sgs)- (-^'^) 



In the first equation, the rate of viscous dissipation of kinetic energy, V-(i>-cr), 
is neglected under the assumption that the flow is virtually unaffected by mi- 
croscopic viscosity on length scales greater than the grid resolution A. This is 
a valid approximation for A ^ t/k, where r/K is the Kolmogorov scale of vis- 
cous dissipation. For the numerical simulation discussed later-on, A ~ 10^ cm, 
whereas r/K ^ 1 cm pi] . The symbolic terms in equation (|17|) account for the 
diffusion, production and dissipation of SGS turbulence energy (see 15 for the 
exact definitions) . Energy transfer from resolved toward subgrid scales is given 
by the rate of production Sggg. The non-local transport term 2) sgs accounts for 
the redistribution of turbulence energy by subgrid scale velocity and pressure 
fluctuations. Furthermore, there are two contributions to the rate of dissipa- 
tion: pe^gs is caused by the viscosity of the fluid, while /oAggs is due to com- 
pression effects. In fact, all of the SGS dynamical terms are non- computable in 
terms of resolved quantities. This is a consequence of the non- linear structure 
of the hydrodynamical equations, which prohibits complete decomposition. In 
consequence, one must find heuristic approximations in terms of computable 
quantities, which are commonly known as SGS closures. 

We apply the customary turbulent-viscosity hypothesis for the rate of pro- 
duction (section 10.1 in (SHI)) the gradient-diffusion hypothesis for turbulent 
transport (section 4.3 in |22]) and the dimensional closure for the rate of vis- 
cous dissipation (section 13.6.3 in [Hni)- Since pressure effects are small for de- 
flagration in degenerate matter, a rather crude closure for Aggs is utilised |4nj . 
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The final result is the fohowing dynamical equation (section 3.1.4 in jl5j): 

,3/2 

r^SgS 



A, 



off 



In particular, an expression analogous to the viscous dissipation term in the 
Navier-Stokes equations, aik = pvSik is substituted for the anisotropic part of 
Tik- The rate of SGS turbulence production is then given by 



^sgs — TikSik — P I 2^sgs|5' I ~ o^sgst^ ) ; (18) 



where fggs = Cyl\eskl^ is the SGS turbulence viscosity. The rate- of- strain 
tensor Sik is the symmetrised spatial derivative of the velocity field: 

The trace of this tensor yields the dilatation of the velocity field, d = Su. 
The scalar 15"*! is formed by total contraction of the trace- free part of the 
rate-of-strain tensor, 



1^1 = VWSfk = y 2 [s^kS^u - \d'^ ■ (20) 

The length scale Aeg is an effective scale of the finite- volume scheme, namely, 
the PPM. The ratio (3 = Aeff/A specifies the smoothing of the fiow on the 
smallest resolved scales due to numerical dissipation. In ^J, we propose a 
method of calculating f3 from numerical realisations of isotropic turbulence. 
For moderately compressible flows, it appears that (5 ~ 1.6. 
Alternatively, a dynamical equation for the turbulence velocity q^^s = 
sgs can be formulated: 



,2 



^ (21) 
The length scales introduced above are defined as follows: 



4 = C«Aefr/^, £^ = aAeff/V2, 4 = 2V2Acfr/C,. (22) 
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The equation for qsgs can be evolved starting with the initial data (7sgs(a;, 0) = 
for a fluid being initially at rest. Moreover, non-integer powers of qsgs do not 
occur, and the functional dependence on Qggs is advantageous for the discreti- 
sation of the diffusion term. However, numerical errors may arise from the non- 
conservative form of equation ()21() . Since we apply the solution predominantly 
to estimate the turbulent flame speed, this caveat is not of much concern. 

At this point, one is left with the problem of determining the closure pa- 
rameters Ck, Cu, and C\. For isotropic turbulence, approximate statistical 
values from analytic theories or numerical data can be found. We adopted the 
constant parameter C\ = —0.2 [32]) and the turbulent diffusion parameter 
Ck = 0.36 was estimated from inertial-subrange properties of flow realisations 
in numerical simulations of forced isotropic turbulence (section 3.2.4 in jl5j). 

A more sophisticated approach is the numerical in situ computation of SGS 
closure parameters from local structural properties of the flow. The underly- 
ing idea is that turbulence in the inertial subrange becomes asymptotically 
self-similar towards smaller length scales. In other words, mostly turbulent 
velocity fluctuations on the smallest numerically resolved length scales deter- 
mine the local energy transfer towards unresolved scales. This idea initiated 
the development of so-called dynamical procedures for the computation of C^. 
The result is a localised closure for the rate of energy transfer Sggg (section 4.3 
in [22]). For the rate of dissipation, eggs, localised closures have been suggested 
as well. Given the computational difficulties and conceptual shortcomings of 
these closures, we decided to apply a statistical method for calculating time- 
dependent mean dissipation parameters in regions containing fuel, flames or 
ash, respectively. In the following section, we will explain the computational 
procedures for Ci, and in detail. A generalisation including dynamical pro- 
cedures for Cfj or Cx as well would be extremely involved and is likely to be 
infeasible in terms of computational costs. Fortunately, the most important 
contributions to SGS dynamics arise from the production and the dissipation 
terms. 



4.2 The semi-localised model 

In order to extract the small scale velocity fluctuations in a simulation, a test 
filter ( )t is applied. This filter smooths the numerically computed flow over a 
characteristic length At = jT^eS: where the factor 7t > 1. It is then possible 
to compute the turbulence stress associated with the intermediate range of 
length scales A ^ Z ^ Ax: 




(23) 
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Applying the eddy- viscosity closure to the trace-free part of TT(wt, W/t), we have 

rTiv^,Vk)sll^ = pTaAT4/'|5*m|2^ (24) 
where pT = {p)t is the test-filtered mass density, and the rate of strain at the 

fTl 1/2 

test filter level is defined by S-j^ = d(^i{pVj-^)T/ Pt- The expression Cj^Axferp 
has the dimension of viscosity. The kinetic energy A;t is defined analogous 
to equation (|18|) . with the test filter in place of the implicit filter and the 

numerically computed velocity v in place of v. 

Invoking the similarity hypothesis that is equal for the eddy-viscosity 
closure at the test filter level and for the unfiltered SGS turbulence stress, the 
anisotropic part of the rate of production in the localised SGS model is given 
by 



is +h d^ia \S*\' "tK.^.)€' \S*\' J^s 

^2.sgs+ gfcsgsd-^.gsgsl^ I - ^^^^ |^*[T]|2Y^^- ^^^> 

The above closure is basically the result of adapting the Germano-Lilly dy- 
namical procedure for the localised Smagorinsky model to the SGS turbulence 
energy model [131 ■ ™ important difference, however, the eddy viscosity 
closure is applied to TT{vi,Vk) rather than the total turbulence stress at the 
test filter level, TT{v'i,Vk)- The stress tensors are related by the Germano iden- 
tity [SHI: 

rT{Vi,Vk) = {T{Vi,Vk))T + TT{Vi,Vk) (26) 

This modification of the dynamical procedure was proposed by Kim et al. |23j . 
It is supported by results from the evaluation of velocity measurements in 
round jets j44| and was explicitly verified on data from simulations of com- 
pressible turbulence driven by stochastic stirring (section 3.2.2 in 15 ). 

A complication arises from C^, becoming negative in some regions of a turbu- 
lent fiow. This is commonly interpreted as backs cattering, i.e. kinetic energy is 
locally transfered across the cutoff from smaller, unresolved vortices towards 
vortices of size larger than A (section 4.4 in [221 )• Including the contribu- 
tions from backscattering in numerical simulations introduces several difficul- 
ties [ini- Firstly, numerical instabilities might be induced, because backscat- 
tering amounts to negative diffusion. Secondly, the SGS turbulence stresses 
must be coupled to the resolved flow in order to consistently account for the 
conversion of SGS turbulence energy into resolved kinetic energy. This is ex- 
actly what one would do in conventional large eddy simulations. However, in 
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combination with a dissipative finite- volume scheme such as the PPM, includ- 
ing the SGS stress terms in the momentum equation does not do much good. 
Inevitably, the kinetic energy produced by positive SGS stress (corresponding 
to negative SGS viscosity) would be injected into modes corresponding to wave 
numbers near the cutoff. These wave numbers, however, are severely affected 
by numerical dissipation and the fluid motion produced by the inverse 
energy transfer would be rapidly damped out. Consequently, backscattering 
would effectively result in enhanced dissipation, i.e. conversion of subgrid scale 
turbulence energy into internal rather than resolved kinetic energy. For this 
reason, the outcome of suppressing backscattering will be investigated in sec- 
tion 

Finding a dynamical procedure for the parameter of SGS dissipation, C^, 
is yet more demanding. The difficulty of determining stems from the fact 
the rate of dissipation is mostly determined by the fluid dynamics on scales 
much smaller than the numerical resolution. Therefore, a localised similarity 
hypothesis is bound to fail. However, one can invoke a statistical argument. 
Considering a certain region of the flow, the average rate of dissipation in that 
region should be roughly balanced by the mean transfer from larger toward 
smaller scales, if the flow is nearly in statistical equilibrium. Even in developing 
flows, a time-dependent statistical value of can be calculated by means of 
energy conservation. The method is loosely based on the variational approach 
of j46j . where the parameter of dissipation is determined by subtracting 
the test-filtered SGS turbulence energy equation (fT7|) from the corresponding 
equation for the unresolved kinetic energy at the level of the test filter. Rather 
than computing locally, we will determine statistical values evolving in time 
from the spatially averaged energy equations. Upon averaging equation (|17j) . 
one obtains 



The diffusion term cancels out, because integrating the divergence of the dif- 
fusive flux over a domain with periodic BCs yields zero. Furthermore, 



i.e. there is vanishing net advection over the whole domain of the flow. The 
turbulence energy at the characteristic scale of the test filter is defined by 




(27) 




(28) 



=0 




(29) 
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and the corresponding averaged dynamical equation is 

Q / oo oo fTl \ 

-g^{pKsgs + PtKt) = (^TT:{vi,Vk)Slf^' J - {p{Xsgs + Csgs) + Pt(At + ex)) • (30) 



Equations (|27)) and (|29|) in combination with the Germano identity (|26|) imply 
the following conservation law for the mean turbulence energy (-STt) associated 
with the smallest resolved scales: 



^^{Kt) = (tt{v,, vk)si;:^ + {TikhsSi^ - TikS.k) - (pt(At + ex)) • (31) 



Substituting the turbulent-viscosity closures for the various production 
terms on the right-hand side, the above equation becomes 

^(Kt) =^ (/5TaAT7^|5*[Tl|2) - ^ ^Kt^I^^) -(ptAt) + (ptct) 

V ' 

(I) 



(II) (III) 

Analogous to the rate of strain at the test filter level, the divergence 
is given by d^^^ = di{pvi)"[ / p^^. The most significant production term is (I) 
which measures the energy transfer across the test filter scale At- The two 
contributions in term (II), on the other hand, are both related to the energy 
transfer across Agfj, where the first expression is calculated from the test- 
filtered and the second expression from the numerically resolved rate of strain, 
respectively. It appears reasonable to assume that the difference of these two 
expressions is marginal relative to (I) in the case of scaling ratios At / of 
the order unity. Furthermore, the evaluation of (II) is particularly costly due to 
several tensor components which have to be test-filtered. Thus, we neglect term 
(II). For similar reasons and because of the smallness of compressibility effects, 
we drop (III) as well. In conclusion, the rate of dissipation eT is approximately 
given by 



(pTeT) ^ -^(^t) + (/9TaATV^|5*[Tl|2) - ^ (|pT(A:Tt^'^' + At)) • (32) 

The crucial step is to conjecture that the relation between the spatially 
averaged dissipation rate and turbulence energy is similar at the cutoff and 



February 5, 2008 7:7 Combustion Theory and Modelling TurbDefl CTM 



W. Schmidt et al. 17 

the test filter level. This implies 

iPTCT) = ^PT + fcx) - iTpklii^ . (33) 

Note that the expression in parentheses is the total turbulence energy at the 
test filter level. For the pressure-dilatation term At, we set 

\r^^CxkT(P\ (34) 



which is analogous to the closure at the subgrid-scale level. The rate of SGS 
dissipation is therefore given by 




(35) 

As opposed to the statistical values for the SGS parameters for steady isotropic 
turbulence, the above equation yields a spatially constant parameter evolving 
in time. This method of calculating in combination with the dynamical 
procedure for Cu makes up the semi-localised SGSTE model. For the numerical 
implementation, two further modifications were added. 

On account of the anisotropy in the vicinity of a flame front, it seems advis- 
able to average over the principal topological subdomains, namely, the interior, 
the exterior and the interface. The latter is identifled by marking all grid cells 
within a certain maximal distance to those cells in which the level set func- 
tion G swaps its sign. With this procedure, the functions c'f'\t), C^\t) and 
Ce (t) are obtained for the mean dissipation parameters in ash, the burning 
zone and fuel, respectively. In the early stage, burning regions might encom- 
pass only a small volume fraction with relatively high surface to volume ratio. 
Hence, the corresponding spatial averages in equation H35() will not be suffi- 
ciently well behaved at the beginning. Although, the dynamics is dominated 
by the fuel domain at this point, both the enumerator and dominator in equa- 
tion (|35j) are smoothed in time via convolution with an exponential damping 
function in order to remove strong oscillations in the ash and fiame regions. 
The characteristic time scale of smoothing is prescribed by the parameter 
T^. An appropriate choice for the time scale has to be found a posteriori. 
Setting ~ O.IT appears to be a good choice in order to get well behaved 
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functions Ce (t), Ce (t) and Ce (t), without overly damping dynamical vari- 
ations (section 4.3.4 in 

5 Numerical Simulations 

The simplest case study one can think of is the evolution of flame fronts in spa- 
tially homogeneous turbulent flows. To that end, we implemented a stochastic 
driving mechanism for the production of turbulence in a cubic domain subject 
to periodic boundary conditions jl5pi4j . The notion of a stochastic force field is 
outlined below. The computational domain is divided into a lattice of subcubes 
of size L = X/a, where L is the characteristic wavelength of the stochastic 
driving force and X is the domain size. In the centre of each subcube, ther- 
monuclear burning is ignited in small spherical regions located at time t = 0, 
when the stirring force begins to act on the fluid. We chose a = 2, giving eight 
subcubes in the computational domain. This pattern is inflnitely repeated in 
space by virtue of the periodic boundary conditions. Gravity is negligible at 
the scales under consideration (section 2.3.3 in ^Sj). Consequently, there are 
no buoyancy effects and turbulence is only produced by stirring. 

The crucial parameter for the evolution of the burning process is the ratio 
C = siam/y, i-e. the ratio of the laminar burning speed to the characteristic 
velocity of the turbulent flow. Assuming developed turbulence, one can apply 
the Kolmogorov scaling law and estimate the magnitude of turbulent velocity 
fluctuations at a separation of the order to the Gibson length: 



The integral length L and characteristic velocity V specify the largest turbu- 
lent vortices in the flow. Setting v'{Iq) = siamj the scaling law for the Gibson 
length becomes 



Obviously, Iq is very sensitive to value of ^. Possible choices of V are restricted 
by the speed of sound Cg. Both Cg and siam are mostly determined by the mass 
density, and so is the ratio Ig/L. In the following, we will consider two dis- 
tinct cases. For ^ > 0.1 and sufficiently high resolution, the Gibson scale is 
just within the range of numerically resolved length scales. In this case, no 
SGS model is required for the flame dynamics and the propagation speed is 
more or less given by the laminar flame speed. If ^ <C 0.1, on the other hand. 




1/3 



(36) 




(37) 



February 5, 2008 



7:7 



Combustion Theory and Modelling 



TurbDeflCTM 



W. Schmidt et al. 



19 



it is impossible to resolve the flame completely. Then the burning process will 
enter a turbulent regime, in which the flame propagation speed is asymptot- 
ically given the SGS turbulence velocity Qsgs- Prior to the discussion of the 
numerical simulations, we give a brief description of the stirring mechanism 
for the production of isotropic turbulent flow. 

5.1 Stochastic forcing 

The specific driving force f{x,t) is composed in spectral space, using a three- 
dimensional generalisation of the scalar Ornstein-Uhlenbeck process, as pro- 
posed in 14_. The evolution of the Fourier transform f{k,t) is given by the 
following Langevin-type stochastic differential equation: 



df{k,t) = -f{k,t)- + FoJ2[^^] S{k-kji^)P^{k).dWt, (38) 



The second term on the right hand side accounts for a random diffusion pro- 
cess, which is constructed from a three-component Wiener process "Wf The 
distribution of each component is normal with zero mean and variance dt. The 
wave vectors kji^ are dual to the position vectors of the cells in the numerical 
discretisation of the fundamental domain. The symmetric tensor P^(fc) is de- 
fined by the linear combination of the projection operators perpendicular and 
parallel to the wave vector. The components of P({k) can be expressed as 



where the spectral weight C determines whether the resulting force field in 
physical space is purely solenoidal, dilatational or a combination of both. The 
variance cr'^{k) specifies the spectrum of the force field. We use a quadratic 
function, which confines the modes of the force to a narrow interval of 
wavenumbers, k £ [0, 2A;o] . The wave number kQ determines the integral length 
scale of the flow, L = 27r/ko. 

The root mean square of the specific driving force is determined by the 
characteristic magnitude Fq and the weight C: 



jlm 




(P,,)c(fc) = CPijik) + (1 - C)Pl^-(fc) = CSij + (1 - 20 




(39) 



/rms = Y^^fjlmit) ' f ,Ut)) ^ {I - 2C + K^)FI 



(40) 



jlm 



Since Fq has the physical dimension of acceleration, it can be expressed as 
the characteristic velocity V of the flow divided by the integral time scale T, 
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which is given by the auto-correlation time T of the driving force (|38|) . Setting 
T = L/V, we have Fq = V/T = LV^, and, starting with a homogeneous fluid 
at rest, the flow is developing towards a fully turbulent steady state within 
about two integral time scales. 

5.2 Quasi-laminar burning 

To begin with, we shall consider the case ^ ~ 1. Then the laminar flame 
propagation is fast enough to burn the smallest numerically resolved eddies 
in less than a turn-over time. Consequently, subgrid scale turbulence does 
not affect the flame dynamics. An estimate of the characteristic velocity V 
for given numerical resolution and laminar burning speed is readily obtained 
from relation ()37() . The effective range of length scales which can be resolved is 
roughly given by L/Aefj = N/aP, where N = A/X is the number of numerical 
cells in one dimension. For Iq ~ A^s, we therfore must have 



For the simulation, which will be discussed in the following, we used 432^ 
grid cells. Setting = 432, relation H41() implies that at most V ~ 4siam is 
admissible. Given a moderate mass density, this would entail an extremely 
low Mach number. However, computing an almost incompressible flow with 
the PPM would be infeasible. On the other hand, for an initial density po ~ 
2.90 • 10^ g cm~^, one obtains siam ~ 1.05 • 10^ cm s~^ through interpolation of 
numerical data taken from Timmes and Woosley [S]. The speed of sound for 
this density is cq ~ 9.70 • 10^ cms~^. Choosing V = 4siam ~ 4.20 • 10^ cms~^ , 
the characteristic Mach number is V/cq ~ 0.043. This is quite small, but 
still computationally manageable with a fully compressible hydro code. The 
resulting Gibson length, Iq ~ 3.3A, allows for some margin between Iq and 



Actually, Landau-Darrieus instabilities would induce a small-scale cellular 
flame structure j47j . Due to the significant numerical dissipation at length 
scales / < lOA jlj, the cellular structure will inevitably get smeared out if 
/q ~ A. For this reason, an effective cellular propagation speed Sceii slightly 
larger than siam would be the correct intrinsic propagation speed in place of 
the laminar burning speed 23|- However, this effect is ignored, as the change 
of the Gibson scale due to the difference between siam and Sceii is only about 
a factor of two. Consequently, using Sceii as intrinsic propagation speed would 
not change the flame dynamics dramatically. 

Regarding the numerical distortion introduced by the passive implementa- 




(41) 



Aeff 



1.6A. 
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tion, one should be on the safe side for mass densities larger than lO^gcm"^. 
Apart from that, the randomisation caused by turbulence tends to diffuse any 
numerical artifacts. In this respect, propagating symmetric flame fronts, say, 
nearly planar or spherical flames, is a more demanding task. Furthermore, 
flow maps prepared form the simulation data clearly show a tight correla- 
tion between the shape of the front and the flow structure (see figure IHI). If 
there were significant spurious propagation or deformation, the evolution of 
the front should become increasingly uncorrelated to the flow. In conclusion, 
the simulations which will be discussed subsequently are likely to give a sound 
description of the flame dynamics, albeit the shortcomings of the level set 
method in the passive implementation. 

The progression of the deflagration in the course of the simulation is illus- 
trated by the sequences of contour plots for the specific internal energy, the 
mass density and the rate-of-strain scalar in figures |2l and El respectively. 
The contours of the zero level set are visible as thin white lines which sepa- 
rate dark regions containing unburned matter of low energy density from the 
brightly coloured regions containing processed material of high energy density. 
In the course of the first integral time T, the regions of burned material are 
expanding gradually. At the same time, they are stretched an folded by the 
solenoidal large-scale flow. In the following second integral time, vortices are 
generated on small scales. Subsequent to the peak of (B) at i = t/T 1.35, 
ash encloses fuel rather than the other way around, and the density of the en- 
closed fuel is noticeably larger than the average density. Around t ~ 1.5 most 
of the fuel has already been consumed by the burning process, and the last 
fuel patches are disappearing quickly. Thus, the peak of burning is reached 
before turbulence is fully developed and the front propagation is affected only 
little by small-scale velocity fluctuations in this simulation. Consequently, we 
refer to this mode of burning as quasi-laminar. 

In panel (b) of figure |3 the mean burning rate {B) is plotted on a dimen- 
sionless scale. One can see that the rate of burning increases exponentially 
in the interval 0.3 ^ t < 1.2. The norm IS"*! of the trace-free rate of strain 
defined in equation (|2()|) is a so-called structural invariant of the flow. Contour 
plots of 1 5* I are shown in figure |31 Regions which are subject to intense strain 
correspond to steep velocity gradients and appear bright in the contour plots. 
These regions tend to form vortices and clearly influence the morphology of 
the flame fronts. As one can see in figure 01 the white lines indicating the zero 
level set tend to be aligned with structures associated with large strain. The 
statistics of \S*\ and two further structural invariants, namely, the vorticity 
uj = |V X i)| and the divergence d = V • i), is plotted in panel (d) of fig- 
ure |3 Because of the small characteristic Mach number of the flow, we have 
d <C a; ~ IS"*!, where the equality of vorticity and rate-of-strain scalar holds 
asymptotically in the limit of incompressible flow. The graphs show that the 
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dns_energy_evol . png 

Figure 1. Three-dimensional simulation of thermonuclear deflagration in a cube with the flame 
propagation speed being equal to the laminar burning speed. The initial density of the C+O fuel is 
po ~ 2.90 ■ lO" gcm^'\ and V = 4siam (C = 0.25). The characteristic Macli number of the fully 
developed turbulent flow is V/co 0.043, where co is the initial sound speed. Shown are 
two-dimensional contour sections of the normalised specific energy e = b/cq at different stages of 

the burning process. 



root mean square (RMS) of 1 5**1, grows exponentially from the first few tenths 
of an integral time scale up to t/T ~ 2, where the stagnation of the growth 
marks fully developed turbulence. A comparison between panels (b) and (d) 
suggests that the growth of the burning rate prior to the maximum correlates 
with the exponentially increasing (IS"*]^)^/^. This underlines the above state- 
ment about the influence of strain onto the flame evolution. The corresponding 
evolution of the RMS momentum and Mach number is plotted in panel (a). 

The rate of change of the mean mass fraction of fuel, {X{C + O)), is also a 
measure of the burning speed. For an energy release e^uc per unit mass, the 
total nuclear energy generated in the whole cubic domain per unit time can 
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Figure 2. 2D contour sections of the relative density fluctuations (p — po)/po corresponding to the 

panels shown in flgurelTI 



be expressed as 

{aLf{B) = -enucM,^h{X{C + O)), (42) 

where M^uh = Poicil-')^ is the total mass contained in the computational do- 
main. On the other hand, the rate of energy production is related to the total 
surface area of the flames, Ap, and the laminar propagation speed, provided 
that compression effects are neglected: 

{aLf{B) = Poenuc^FSlam- (43) 



Combining equations and with M^ub = 8poL^, the approximate total 
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dns_strain_evol . png 

Figure 3. 2D contour sections of the logarithmic dimensionlcss rate of strain, log]^Q{T|S'*|), 
corresponding to the panels shown in figures [Tl and 1^ 



surface area is given by 



Ap~l^(X(C + 0)). (44) 

■^lam 



The graph of the normahsed surface area, 



■«^=8^ = -^^£(^-*(C + 0)), (45) 



is shown in the panel (c) of figure HI The exponential growth of the burning 
rate is manifest in this plot as well. At the peak, Ap ~ 1, which verifies that 
the flames experience only little wrinkling due to small vortices. This result 
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Figure 4. Evolution of dimensionless statistical moments for the simulation illustrated in the 
series of figures [Tl 1^ and ISl The top panels show plots of the RMS force, momentum and Mach 
number as well as the average nuclear energy generation rate in combination with the chemical 
composition. Furthermore, a measure for the mean flame speed and the averaged structural 
invariants of the flow are plotted in the panels at the bottom. 

agrees with the impression of rather smooth flames in figure ^ Also plotted is 
the graph of —{TX{C + 0))/{X{C + O)), which is a measure of the ratio of 
the flame surface area to the amount of still unburned material.^ The mean of 
this ratio follows a nearly exponential law even beyond the peak of {B) and, 
thus, can be considered as an invariant measure for the burning intensity even 
when the flames are already diminishing. 



^Strictly, the volume of fuel left at a certain time would be given by {aLY^ {pX(C + O)) / po ■ However, 
the mass-weighted fraction of C+O was not calculated in the simulation. 
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5.3 Burning dominated by turbulence 

In the case ^ <C 1, i.e. the laminar burning speed is smah compared to the 
chacteristic velocity of the flow, the range of length scales between the Gib- 
son scale and the integral length scale becomes very large. Consequently, it 
is impossible to resolve the flame dynamics completely. For example, setting 
po ~ 2.90 • 10^ gcm~'^, which is by an order of a magnitude smaller than the 
density chosen in section EH yields siam ~ 9.78 • 10^ cms~^. Choosing a char- 
acteristic velocity V = lOOsiami we have the Mach number V/cq « 0.15, and 
the Gibson scale becomes ~ 10~^L. Obviously, Iq <^ A for any feasible 
numerical resolution. A subgrid scale model is therefore mandatory. In this 
section, several simulations of thermonuclear deflagration in the cube with the 
turbulent flame speed given by equations and ©, respectively, are dis- 
cussed. The SGS turbulence velocity q'sgs is computed via equation (PT|) . The 
initial mass density and the characteristic velocity are as specified above. Oth- 
erwise, the same parameters as in section are used. However, the resolution 
is reduced by a factor two. So there are 216^ grid cells. 

A couple of simulations were performed with the semi-localised model and 
the maximum relation (jHJ with Ct = 1. In one case, we coupled the SGS 
stresses to the resolved flow and included backscattering, whereas no coupling 
was applied and backscattering was suppressed by setting the SGS viscosity pa- 
rameter equal to Cj" = max(0, C^) in the other case. The validity of neglecting 
the SGS stress terms in the momentum equation (|14|1 has been investigated in 
several hydrodynamical simulations with PPM |48|l41j. For combustion prob- 
lems, the SGS model then runs in a passive mode and provides the turbulent 
flame speed only. 

If negative values of C^, are admissible, however, the SGS stress terms must 
be included, because otherwise backscattering would convert SGS turbulence 
energy into heat rather than kinetic energy on resolved scales. Hence, backscat- 
tering necessitates an active SGS model. For a fully consistent treatment, the 
terms v • (V • r) and pcsgs have to be added on the right hand side of the 
conservation law for the total energy etot = + ^int) which account for the 
transfer of kinetic energy between resolved and subgrid scales and the pro- 
duction of internal energy due to the viscous dissipation of SGS turbulence 
energy, respectively. 

On the other hand, if backscattering is suppressed, a considerably simplified 
scheme is applied, where the dissipation of kinetic energy is solely of numerical 
origin, and Qsgs is treated as a passive scalar. In order to account for the 
exchange of energy between the resolved total energy, etot = I'W^ + Cint, and 
the SGS turbulence energy in a rudimentary fashion, the Lagrangian rate of 
change of /cggs is subtracted from the the conservation law for the total energy. 
Locally, this introduces a certain error due to the diffusive transport of SGS 
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turbulence energy. In fact, the transport term in equation (|18|) changes the 
unresolved energy without affecting the resolved energy budget. We ignore this 
contribution in the energy update, if backscattering is suppressed and the SGS 
model is not completely coupled to the resolved flow. In the case of complete 
coupling, however, the exact SGS transfer and dissipation rate are accounted 
for in the dynamical equation for eint- Of course, complete coupling would seem 
appropriate regardless of the treatment of inverse energy transfer. However, if 
backscattering is suppressed, we found that abandoning the SGS stresses and 
using the approximate energy update as outlined above changes the results 
only little, while increasing the computational efficiency significantly. 

The test filter for the computation of the production parameter in the semi- 
localised SGS model is numerically implemented as follows. Orthogonal one- 
dimensional mesh filters with nine supporting nodes are applied in the direc- 
tion of each coordinate axis. The filter weights are determined by matching 
the Fourier transform of the kernel as accurately as possible to the spectral 
representation of an analytic box filter. Then a single free parameters remains, 
which is the filter scaling ratio jt- For a given number of supporting nodes, 
A'^x, an optimal value of 7t can be found from further constraining the Fourier 
transforms of the mesh and the analytic box filter, respectively, to be equal at 
the cutoff wavelength (appendix A. 1.1 in 15 ). In the case iVx = 9, we have 
At 3.75Aeff « 6.74A. Although 7t ~ 3.75 is considerably larger than a fac- 
tor of two, which is commonly suggested in the literature, we obtained optimal 
results with this setting rather than with test filters of smaller characteristic 
length (section 4.3.3 in |T5]). 

Statistical results from the simulations are shown in figures El and El The 
evolution of the mean burning rate {B) and the corresponding fuel consump- 
tion together with the helium and nickel production is shown in the bottom 
panels (g,h,i) of figure [3 For comparison, the RMS forcing, momentum and 
Mach number are plotted in top row of panels (a,b,c) for each simulation. In 
the course of the first integral time, the slope of {B) in logarithmic scaling is 
rather slowly rising. This indicates predominantly laminar burning. The oscil- 
lations of {B) at early time are caused by numerical discretisation errors, since 
the burning regions initially tend to become elongated into thin shapes and, 
in consequence, are only marginally resolved. As the flow becomes increas- 
ingly turbulent and a growing fraction of the total flame surface is subject to 
an enhanced propagation speed st ^ sianu the rate of burning rises rapidly. 
Eventually, the phase of exponentially growing energy release passes over into 
fading combustion once the greater part of the fuel has been exhausted. The 
transition point between the quasi-laminar and the turbulent burning phase 
can be estimated from the tangents to the almost linear portions of the graph 
of {B) in logarithmic scaling. By means of the plot of the mass-weighted sta- 
tistical moments of Qsgs in the panels (d,e,f), one can see that the transition 
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Figure 5. Evolution of dimensionless statistical moments for simulations with different flame 
speed models and SGS closures. The normalised RMS of the specific stirring force and the resolved 
momentum of the flow are shown in the top panels. The first three moments of the SGS turbulence 
velocity are plotted in the middle row of panels. Note that the mass-weighted mean of qegs is scaled 
in units of the laminar burning speed siam- In the bottom panels, the mean burning rate and the 
average mass fractions of fuel and processed nuclei are shown. 



coincides with (/5(?sgs)/(poSiam) ~ 3. The somewhat larger threshold value of 
the mean SGS turbulence velocity relative to the laminar burning speed for 
the onset of rapid turbulent burning is possibly a consequence of intermit- 
tency. This is also indicated by the large standard deviation, a{pqsgs), which 
is comparable to the mean, (pqsgs), during the production phase. The mass- 
weighted skewness, skew(/9(7sgs), is particularly large in the early phase when 
eddies are forming locally, whereas it approaches a value near unity in the 
regime of statistically stationary and homogeneous turbulence. 

Comparing the left and the middle column of plots in figure |SJ the evolution 
of the burning process in the simulations which differ only by the coupling of 
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Figure 6. Comparison between two different variants of the semi-localised SGS model. 
Backscattering is included in one case, whereas it is suppressed in the other case. The top panels 
show plots of the mean production parameter C'v in regions containing ash, flames and fuel, 
respectively. The middle panels show the corresponding plots of the dissipation parameter, and the 
different contributions in equation I21i for the time evolution of gsgs are plotted in the bottom 

panels. 



the SGS model to the momentum equation and the treatment of inverse energy 
transfer appear quite similar. This can be seen in figures [7| and as well, which 
show contour plots of the total energy per unit mass at different instants of 
time for both simulations. However, the burning process proceeds faster in the 
simulation without backscattering in comparison to the simulation with the 
fully coupled SGS model. Correspondingly, the peak of (B) is delayed in the 
latter case. 

The differences in the SGS model are illustrated in figure IHl The averages of 
the production parameter Ci, in regions containing, respectively, ash, flames 
and fuel are plotted in panels (a) and (b), respectively. If backscattering is 
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Figure 7. Simulation of thermonuclear deflagration for an initial density po 2.90 ■ 10* gcm~^. 
W*lam = 100, corresponding to a characteristic Mach number V/cq ^ 0.15. The turbulent flame 

speed is given by st = max(siami 9sgs). The turbulent velocity gsgs is computed with the 
semi-localised SGS model, where = max(0, C^) is set as parameter of turbulence production. 
Shown are 2D contour sections of the normalised specific energy e = e/cg at different stages of the 

burning process. 



included, the mean of Ci, within the fuel is initially small. In the course of 
turbulence production, the parameter is growing and, eventually, (Cu) ^ 0.04 
in the statistically stationary regime. On the other hand, if backscattering is 
suppressed, the mean of Ci, changes only little in the dominating topological 
region. The local values of d, are fluctuating considerably for t < 1. This can 
be seen from the strongly oscillating averages of Ci, in ash and flames, which 
initially fill quite narrow spatial regions and encounter varying conditions while 
being advected by the flow. The time evolution of the dissipation parameters 
shown in panels (c) and (d), on the other hand, exhibits more or less the same 
trend. During the first turn-over time, Ce vanishes identically. At time i ~ 1, 
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Figure 8. 2D contour sections of SGS turbulence velocity relative to the laminar burning speed, 
gsgs/sjami in logarithmic scaling. The panels correspond to those shown in Figure 171 



small turbulent vortices begin to form and turbulent dissipation sets in. In 
statistical equilibrium, assumes a nearly constant value of about 0.65 for 
the fully coupled model and 0.75 without backscattering. In the latter case, 
a larger dissipation rate compensates the suppressed inverse energy transfer, 
as one can see from the plots of the mean rate of production and dissipation 
corresponding to the source terms on the right hand side of equation H21|) in 
panels (e) and (f) of figure IHl It is obvious that complete coupling significantly 
reduces the rate of turbulence production. Nevertheless, the mean value of 
gsgs is found to be nearly the same in the statistically stationary regime for 
both simulations. As mentioned in section 14. 2| we suspect that turbulence 
production is systematically underestimated by the fully coupled SGS model 
in combination with the PPM, because the kinetic energy injected through 
backscattering into modes of high wave number will be quickly dissipated by 
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Figure 9. Simulation with the same parameter as specified in the caption of figure|7] however, 
with Cjy as production parameter and complete coupling of SGS stresses to the flow. 



numerical viscosity. An impression of the spatiotemporal evolution of ^sgs is 
given by the contour plots in figures |H1 and ^1 respectively. 

On account of results from laboratory measurements, Kim et al. argue that 
Pocheau's relation Q with n = 2 and Ct = 20/3 gives the most accurate pre- 
diction of the turbulent flame speed The outcome of running a simulation 
with the flame speed relation 



St = SlamJl + f f^V, (46) 
V -3 V'Slam/ 



is demonstrated by the statistics in the colum of panels (c,f,i) on the very right 
of figure El Backscattering is also suppressed in this simulation. Now the peak 
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Figure 10. 2D contour sections of SGS turbulence velocity relative to the laminar burning speed, 
gsgs/siam in logarithmic scaling. The panels correspond to those shown in Figure 1^ 



of the burning rate is reached even faster than in the case of the simulation 
with Ct = 1. In fact, the bulk of the burning takes place when the level of 
SGS turbulence is quite low. Accordingly, the plots of contour sections of the 
specific energy in figure ^2 show that the flame surface is smoother and less 
corrugated by the flow in the course of the burning process. The slope of {B) is 
steepening significantly just for {pqsgs) / {poSia,m) ~ 1- This would suggest that 
Ct = 20/3 is, indeed, a feasible choice. On the other hand, Ct = 4/3 is favoured 
by Peters which yields more or less the same behaviour as in the case 
Ct = 1. From our current understanding, we should consider Ct as a parameter 
of the flame speed model, which is to be chosen within reasonable limits and 
validated a posteriori by the results obtained for a particular application. 
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Figure 11. The physical parameters are the same as for the simulations in the previous figures, but 
this time the flame speed relation is given by st = siam[l + ^(9sgs/siam)^]^ ■ Shown are 2D 
contour sections of the normalised specific energy e = e/cg at different stages of the burning 

process. 



6 Conclusion 

The numerical simulation of thermonuclear deflagration in a box subject to 
stochastic stirring was utilised as a test problem for the study of flame speed 
models. The evolution of the flame front was computed by means of the level 
set method in the so-called passive implementation. Essentially, an effective 
flame propagation speed must be calculated, if the Gibson scale is small com- 
pared to the resolution of the computational grid. A subgrid scale (SGS) model 
based on the budget of turbulence energy determines a velocity scale which 
is proportional to the propagation speed of flame fronts in the fully turbulent 
regime. Some of the closure parameters of the SGS model are locally calculated 
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with dynamical procedures. Thus, we have a semi-locahsed model. 

Particularly, we compared two variants of this model. In one case, inverse 
energy transfer from subgrid toward resolved scales was included, in the other 
case it was suppressed. Inverse energy transfer is also known as backscattcring. 
For a consistent treatment of backscattering, complete coupling of the SGS 
model and the resolved hydrodynamics is indispensable. In combination with 
the piece-wise parabolic method, this entails difficulties stemming from the 
significant numerical viscosity of the scheme. But we obtained sensible results 
when suppressing backscattering and applying a simplified SGS model with 
partial coupling. 



Depending on the constant of proportionality \/Ct in the asymptotic flame 

speed relation, we found the transition from laminar to turbulent burning at 
noticeably different points in the course of turbulence production. This tran- 
sition comes about once the SGS turbulence velocity must exceed the laminar 
burning speed in a significant volume fraction of the computational domain. 
Then the nuclear energy generation grows at a much higher rate, and the 
flame surface develops an intricate structure due to the stretching and folding 
caused by turbulent vortices. If Ct is about unity, the peak of nuclear energy 
release appears roughly when the turbulent flow becomes statistically station- 
ary and homogeneous. For larger values of Ct, most of the fuel is consumed in 
advance of turbulence becoming fully developed. We propose to consider Ct 
as a control parameter, which regulates the overall rapidness of the burning 
process. 

The semi-localised SGS model presented here is especially suitable for any 
kind of transient and inhomogeneous turbulent combustion process. It is the 
first implementation of this kind of SGS model for an astrophysical application, 
namely, the numerical simulation of thermonuclear supernovae. 
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